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Abstract We compare the performance of two alternative algorithms which aim to construct a force- 
free magnetic field given suitable boundary conditions. For this comparison, we have implemented both 
algorithms on the same finite element grid which uses Whitney forms to describe the fields within the 
grid cells. The additional use of conjugate gradient and multigrid iterations result in quite effective 
codes. 

The Grad-Rubin and Wheatland-Sturrock-Roumeliotis algorithms both perform well for the recon- 
struction of a known analytic force-free field. For more arbitrary boundary conditions the Wheatland- 
Sturrock-Roumeliotis approach has some difficulties because it requires overdetermined boundary infor- 
mation which may include inconsistencies. The Grad-Rubin code on the other hand loses convergence 
for strong current densities. For the example we have investigated, however, the maximum possible 
current density seems to be not far from the limit beyond which a force free field cannot exist anymore 
for a given normal magnetic field intensity on the boundary. 

Keywords Coronal magnetic field, force-free field extrapolation 



1 Introduction 

With the advent of vector magnetographs which measure the line-of-sight component of the photo- 
spheric magnetic field and, except for a 180° ambiguity, also its component normal to the line-of-sight, 
the interest in extrapolating these measurements into the corona have grown enormously. 

The line-of-sight component of the photospheric magnetic field has been observed for decades now 
but these measurements alone supply only boundary information at best sufficient for a Laplace field 
model of the coronal magnetic field. The vector magnetograph observations now available considerably 
constrain the photospheric horizontal field and allow estimates also of the coronal current density. As 
a consequence much more realistic coronal field models can be based on these observations. 

Since the magnetic field, at least in the lower corona, completely dominates the plasma forces, a 
"force-free" approximation of the field for stationary situations seems to be a tolerable assumption. 
"Force" in this context means the magnetic Lorentz force j xB. All other MHD forces like gravity, 
pressure, etc. are neglected because they are about three orders of magnitude smaller than jB. Hence, 
the current and magnetic field vectors should be aligned to better than half a degree. 

These simplifications accepted, the magnetic field in some domain V of the corona may be described 

by 

V-B = 0; VxB^j; jxB==0. (1) 

The alignment of current and magnetic field causes the problem to be nonlinear, hence the questions 
which boundary information is to be supplied and how to solve ([T]) are by no means trivial. 
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Boundary conditions which seem necessary and sufficient for ([T]) are ( Boulmezaoud and Amari] 

l2nonh 

n-B on dV : n-VxB on (dV)" or {dV)+ , (2) 

where {dV)^ is that part of the surface of V where n-B is either > or < 0. But not all the boundary 
values which comply with ([2]) are allowed. The field line twist which can be stationarily maintained, 
i.e., n-VxB on dV, is limited by the field energy sustained from the exterior, i.e., by n-B (see section 



Quite some effort has gone into attempts to solve ([T]) for given boundary values either in the form 
^ or differently. Among the most promising schemes are those suggested a long time ago by Grad 
and Rubin l(l958i . abbreviated GR) and more recently by Wheatland, Stur rock, and Roumeliotis ( 2000l . 



abbreviated WSR ). The GR code has firs t been implemented by Sakura i ( 198lh and further developed 



by Amari et al. ((l999t ). R egnier et al. (j2002l ) and Wheatland (2004'). The WSR scheme has been 



extended by Wiegelmann ( 20041 ) and Wiegelmann and Inhester (|2003f ). The development still is an 



active area of research. A compa rison and description of codes in use and other alternative approaches 
can be found in Schrijver et al. (I2006D . 



A difficulty for a thorough comparison of the various schemes, however, is the fact that they are 
often implemented very heterogeneously and also include many different details which could easily 
speed up or slow down their performance and sometimes may even obscure the basic advantages or 
disadvantages of an approach. 

In order the make the two schemes comparable, we apply them to the same problem defined on 
exactly identical grids. Error norms to measure the performance are exactly the same. The special grid 
which we use is explained in section [2l We are convinced that it has many advantages for electromag- 
netic problems like ([T]). 

The WSR and the GR algorithms have been extensively described elsewhere, so that in sections |3] 
and[4]we restrict ourselves to the basics and rather emphasize some or the details of our implementation. 
The numerical results obtained for two different examples are presented and discussed in the final parts 
[5]and[i 



2 The grid 

For the problem we want to investigate the choice of the grid and the representation of the fields is 
very crucial. We found most suitable for our purposes a finite element grid which allows to transform 
standard vector analysis consistently into discrete space. It is relat ed to finite difference grids with 
staggered field components and similar grids based on Yee's scheme ( Yee .11 19661 ) . In part of the math- 



ematical literature these special finite elements are called discrete Whitney forms (Bossavit, 19881) 
because they have very much in common with continuous differential forms. In fact, some of the finite 
elements have been known for a long time, however, the way they are related among each other by 
differentiation operations and to their dual sp ace analogues is relatively new and a matter of cur- 
rent research (iHiptmair .1120011 : lGradi naru."2002^ . The elements are particularly suited for a numerical 



treatment of electromagnetic fields ( Teixeira.i.200l] ) 



We here use the elements in the most simple form to lowest order and on a regular cubic grid which 
spans our computational domain, a square box V = [0, 1] x [0, 1] x [0, 1]. The nxnxn grid cells each 
have a size h — 1/n. The cell vertices are located at ih, i = 0,n; the cell centres are at (z — 0.5)h, 
i = l,n. 

In Figure [1] we present the four types of finite elements which we use here. Each form a Hilbert 
space of either scalar or vector valued piecewise first or zero order polynomials determined by a set 
of parameters which are representative averages of the field to be described. The respective averaging 
areas are indicated in colour. 

Lagrangian elements (0-forms) have as parameters the function values at the 8 vertices of a cell. 
The element function is linear insid e the grid cell s uch that the correct values are met at the vertices. 

The Nedelec elements fl-forms. lNedelec.1ll986[ ) are a discrete representation of a vector field. Its 
element parameters are the averages of a field component over a cell edge along the respective com- 
ponent direction. The element function is constant along the edge and varies linearly in transverse 
direction matching the right values at the four cell edges of the given direction. 
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0-form l-form 2-form 3-form 

Lagrange element Nedelec element Raviart- Thomas ele- FV element 

for scalars £ R for vector fields ment for vector fields for scalars G R 



Fig. 1 Lowest order Whitney forms for a square grid cell. 



primary grid box 




Fig. 2 Dual grid cell in relation to a primary (coloured) grid cell. 

The Raviart-Thomas element (2-forins. lRaviart and Thomas.|[T977l ) has as element parameters the 
face averages of the field component normal to the face. The element function varies linearly in this 
normal direction and is constant across the face plane. 

The last element we need is a finite volume element (3-form) for a scalar function approximation. 
It has only a single parameter per grid cell which represents the average of the scalar over the entire 
cell. 

Every vector differential operation transforms an n-form element in a natural way into a n+ l-form 
element: 

0-form — >- l-form 2-form — div^ 3-form 



As for continuous differential forms, a double differentiation gives exactly a zero field and insures thus 
that curlograd and divocurl vanish identically also for the discrete forms. This vanishing of double 
differentiation is not a consequence of the precision with which we approximate the differentiation of 
discrete forms but is due to th e fact that th e boundary of a boundary is an empty set, hence is a 
consequence of geometry alone ( Janich.ll200lh . It is therefore not surprising that this rule holds also 
for discrete forms. 

In order to allow for Laplacians, the alternative combination of vector differential operators, we 
need to introduce the dual grid as opposed to the primary grid described above. The dual grid is in 
case of the regular square primary grid again a square grid but shifted by half a grid size in each 
axis direction such that vertices of the dual grid are located at the cell centres of the primary grid 
(see Figure [2]). Hence there is a natural association between primary 3- forms and dual 0- forms. This 
association also holds for the other forms, between the n- forms of the dual grid and primary (3 — n) - 
forms. For general grids this association is manifested by the Hodge (or ★) transform ( Hiptmair]|200ll ). 
For the square grid we use here, the Hodge transform is simple because the finite element parameters 
of the dual form are the same as the corresponding parameters of the primary element (except for 
domain boundary effects). Note, however, that while the finite element parameters remain unaltered 
under this transformation, their interpretation and their functional representation changes. 

Of course, the forms of the dual grid are connected among each other by differentiations in the 
same way as for the primary grid. The final pattern of forms together with the mappings among them 
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Fig. 3 Finite element parameters entering into the calculation of the z component of j x B at the cell centre. 
The red arrows on the cell edges denote magnetic field components, the blue face centered arrows the current 
components involved. 
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In this scheme, the usual 6-point stencil of a discrete Laplacian operated on a O-form can be realized 
by divTkrgrad where the star denotes the Hodge transform and the result is a dual 3-form. Likewise, the 
Laplacian on a 1-form can be written as div^grad — curl*curl and returns a dual 2-form. Except for 
boundary effects, the primary and dual grids are on equal footing. 

For our problem, the magnetic field B is considered a primary 1-form (or dual 2-form) . The current 
density j then is a primary 2-form. With these prescriptions, we can perform all differential operation 
on scalar and on vector fields in a consistent way. 

We just mention in passing that these forms can be equally well constructed on an irregular grid 
and element functions can be higher order polynomials if more element parameters are adequately 
provided. Hence high order difference formulas can be set up in a systematic way. Multigrid extensions 
of the Whitney forms are a matter of current research (|Gradinaru.ii2002l) . 



3 The Wheatland-Sturrock-Roumeliotis (WSR) algorithm 



The algorithm proposed by Wheatland, Sturrock, and Roumeliotis ( 2000[ ) tries to find the force-free 



field B from the argument which minimizes a penalty function L, i.e., 

B = argmin(L), L{B) ^ J \wjxB\^ + J |V-Bp. (3) 

V V 

In fact, the integrals can be looked upon as a Hilbert product on the respective finite element space 
defined in the previous section. This conception helps greatly when programming L and its derivatives. 
The calculation of V • B and V x B throughout V requires the knowledge of n • B and n x B on the 
whole of the surface dV. This is more than ([2]) prescribes and the problem of minimising L under 
these restrictive conditions is clearly overdetermined. If inconsistent boundary conditions arc imposed 
on a decent minimum may never be reached. We therefore allow for the option in our program 
to vary the normal and/or tangential field components on individual faces of V. This essentially is 
equivalent to setting n ■ B from V • B = and n x B consistent with a vanishing Lorentz force on the 
respective boundary face. 

In this sense the second integral in ^ is calculated from the squared (dual) 3-form V-B defined 
at the cell vertices by summing over all vertices. Vertices on a domain surface, edge or corner are 
especially weighted with factor 0.5, 0.25 and 0.125 respectively. 

The j x B integral is a little more involved. Formally the exterior product of a 1-form and a 2- 
form should give a 3-form. Hence the x product here has not the property of an exterior product in 
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differential geometry. Yet a 3-forni, however one for each component of j x B individually, gives the 
most compact stencil for this expression. In Figure[3]we show the element parameters which are needed 
for the z component of the Lorentz force at a cell centre. E.g., the contribution jx^y is obtained from 
an average over the two a;-faces of the cell. For each of these faces, jx^y is calculated by multiplying 
2x at its centre with the average of the two By values from the two y-edges of this face. The integral 
in Q then is a simple sum over all three components of the Lorentz force. For each component, the 
respective squared 3-form elements are summed over all cell centres. 

In order to obtain argmin i(B), WSR proposed a simple Landweber scheme by iteratively advancing 

B('+i) = b(') + s(5B« , ,5B« = -■ ^(B«) , (4) 

c/B 

which guarantees that L decreases at every step provided step size s is small enough. In this scheme, 
improvements of B are strictly along the negative gradient direction and step sizes are only guessed 
and reduced if necessary. 

We have instead implemented an unpreconditioned conjugate gradient iteration which at every 
iteration step performs an exact line search to the minimum of L\s) = L(B + s (5B) along the respective 
search direction. Moreover, it selects an improved search direction JB instead of the gradient as in (IH) . 

No te that in contras t to ex isting implementations of this scheme bv lWheatland. Sturrock. and RoumeliotisJ 
(12000') or Wiegelmann (2004) who programmed the discretization of the analytical derivative (|4]), we 
always calculate the numerically more consistent derivative of the discretized function L. 

Conjugate gradient solvers are optimal for linear problems for which the objective function L 
depends to second order on the components of B. Our problem, however, is nonlinear and ([3]) is of 
fourth order in B through the j x B term. We make use of our formulation of ([3]) on the special grid 
introduced in the previous section in order to calculate all five polynomial coefficients of L'{s) in one go. 
This enables us to perform the exact line search at every step without much effort by a single function 
call. From the new minimum, the new search direction JB^*^^) is chosen so that it is a descent direction 
and also i7-orthogonal to the previous search direction (5B(*''. H here is the local Hessian d^L/dBidBj 
at the line search minimum. Likewise, we can also choose the Hestenes-Stiefel variant which yields a 
new search direction which is i/-orthogonal with respect to some average Hessian H. 

A parameter still to be determined in fS]) is w. Its choice will be discussed later along with the pre- 
sentation of the results. In fact, it will turn out to be favourable to make it a space-dependent w(x). This 
i s ano ther difference wit h res pect to the implementations of IWheatland. Sturrock. and RoumeliotisJ 
(|2000h and Wiegelmann (j2004l ). who took w a function of B. 



4 The Grad-Rubin (GR) algorithm 

Several variations of this algorithm exist. While the previous approach to find a force-free field reduced 
the problem to a formal optimization procedure, the approach by Grad and Rubin (1958) is inspired 
by a quasi-physical relaxation: at any time in the iteration the current j = a*^"-'B'^"^ produces via the 
Biot-Savart law a new field B'^"+-'^^. According to the differences between B("+^) and B^"\ a is then 
redistributed along the new field lines giving rise to a new current and hence a new Biot-Savart field. 
In our code, we distribute a along given field lines by solving 

B(«).Va(")=0 (5) 

for given B*^"^ and boundary values for a on dV. In the next step we correct the field by solving for 
the vector potential ^ A of the field update 

ASA = 6j, (6) 
where V-5A = and Sj = VxB^") - a(")B(") , 

with boundary conditions nx 5A = and (n • V) (n-(5A) = 0. These boundary conditions insure that 
V- (5 A vanishes also on the domain boundaries. In the subsequent field correction 

B(»+i) ^ B(") + Vx 5A, (7) 

the normal components of B on the domain boundaries remain unchanged. 
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In terms of the forms introduced above, the scheme looks as follows: 

Sj 
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The scheme basically starts form a potential solution on the left and then cycles the square at the centre. 
Any time B has been updated, a is remapped using ^ before the residual current Sj is calculated. 
Due to numerical discretization errors in the integration of a in ([5]), 6j may have a spurious divergence 
which is checked and if necessary, eliminated by iterating V-(5j = a few times while preserving normal 
boundary conditions and V x . 

While the Poisson equation ([6]) can efficiently be solved with a niultigrid solver, the major compu- 
tational effort is spent in solving ([5]) to the required precision. Since ([5|) is first order, boundary values 
need only be supplied at one end of the field line (or characteristic). With ([2]) this is automatically 
satisfied, however this choice of boundary conditions deprives us of the freedom to emphasize those 
boundary areas where we assume the observations to be more reliable. We therefore do not make the 
distinction between {dV)~^ and {dV)~ as in ^ but we safeguard our solver of ([5]) against inconsis- 
tent boundary values for a by attaching a weight with every boundary value for a. The final value 
on the characteristic is the according weighted average from both end points. This way, the influence 
from uncertain boundary values on the side walls or from imprecise measurements on the bottom 
(photosphcric) boundary can be suppressed. 

In principle, (O is solved by mapping every cell centre along a field line to the boundary and 
interpolating the boundary values to its foot point. To minimise the number of field line calculations, 
we store the a value also in every cell the calculated field line intersects along with the intersection 
coordinates. For about 2/3 of the cells, the field line calculation from their centre then can be discarded, 
because they have previously been intersected by so many field lines, that their a value can reliably 
be interpolated form the intersection information stored. 



5 Results 

We have tested our codes with two different model fields. The first is the Low and Lou ( 1990l ) field 
model, one of the few analytic force free field solutions which now has almost become a standard for 
tests of force-free reconstruction schemes. For the second test model, a twisted flux tube, we do not 
have an analytic solution. For the existence of a solution we rely on the symmetry of the boundary 
conditions supplied. 



5.1 Low and Lou model 



As a test case model we use the analytical nonlinear force- free field solution of Low and Lou ( 1990[ ). The 



field results form a multipole with eigenvalue = 0.42659. We placed it at a depth of I — 0.2 below the 
centre of the ground plane and with orientation in the x-z plane with 45 degrees inclination to the x 
and z axes. The field was scaled so that at the ground plane the vertical field strength ranges between 
-9.95 and 5.09, the alpha values between -16.9 and 5.61 and the vertical current density between -67.9 
and 106. A field line plot of this model is shown in the upper left of Figure S) 

The result of the extrapolation is displayed in the other three panels, all for a grid size n — 64. We 
find that the WSR iteration is strongly biased towards regions where the field strength is large. This 
very probably is due to the Lorentz force ter m in (|3l) which increases with O(B^). For this r eason, the 
original object function L proposed by WSR (I Wheatland. Sturrock. and Roumeliotis.ll2000[ ) had w ^ 
1/|B|. This choice, however, makes L a rational function of the B components which we expect to 
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Fig. 4 Representative field lines of the original Low and Lou (1990) model [upper left), the result of the WSR 
iteration with w — const{upper right), the result of the WSR iteration with w = l/|_Bpotl (lower left), the result 
of the GR iteration [lower right). The colour code at the bottom represents Bz, the top plane shows the vertical 
projection of the field lines. 

slow down the convergence of the minimization iteration. We therefore prefer the choice w ~ l/|Bpot| 
where Bpot is the potential field consistent with the normal component of the given boundary magnetic 
field. The potential field is obtained easily and our choice of w similarly reduces the influence of regions 
where the field is expected to be strong. Since w is not changed during the iteration, it does not destroy 
the analytic properties of L(B). 

The effect of this weighting is remarkable, since -Bpot varies by three orders of magnitude between 
0.027 and 28 throughout V. In Figure[4]we show the result of the WSR iteration with w = const (upper 
right) and w = const/ 1 Bpot | (lower left). The discrepancy of the former at larger heights in regions of 
small magnetic field strength is immediately apparent. The lower right panel shows the result of the 
GR iteration. It reproduces the original in weak and strong field regions equally well to at least the 
plotting precision. Details in the original field lines are reproduced exactly. 

For the WSR iteration the question remains how to choose the constant in the weighting factor w. 
In Figure [5] we display the iteration history of the two integral terms in ^ for different constants in 
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Fig. 5 Decrease of the two integral terms in ([3]) during the conjugate gradient iteration for a grid n — 32 (top) 
and n = 64 {bottom). The weight function was chosen w — const/ |Bpot| with different values for the constant 
as marked at the lower left end of the curves. The tick marks along the curve indicate every 50th (n=32) or 
100th (n=64) iteration step. The circle in each diagram indicates the integral values of the analytic model 
calculated on the respective grid. 

w ~ const/|Bpot| and for two different grid sizes n = 32 and 64. The iterations were started with an 
initial B*^*^^ =0 inside V and the right boundary conditions for the normal and tangential fields on 
dV. Note that for B^") the diver gence and the Lorentz force do not vanish because of these boundary 
values imposed. 

A value const ^ 1 seems optimal which is not surprising because with this choice of w both terms 
are of the same order B^/h. Smaller values of the constant lead to a faster reduction of the divB term, 
larger values yield a bias towards the elimination of the Lorentz force. In Figure [5] we have chosen the 
axes scales for the ordinate and the abscissa equally, so that the curve (const ^ 20) which decreases 
closest to 45 degrees eliminates both integral terms at about equal rates. 

For the case w = const similar plots can be produced as in Figure [5] and the optimal constant then 
is around 0.05. It depends, however, on the magnetic field strength of the model since now the two 
terms in (jS]) have different units. 

The number of iterations performed in Figure [5] was 10 xn where n is the number of cells along 
each coordinate axis. From the similarity of both diagrams we conclude that roughly the number of 
iterations needed to reduce L by a given factor increases proportional to n. This result seems to differ 



9 




log<ldivB|-> 



Fig. 6 Same as lower diagram in Figure [5] but with a potential field as initial B for the iteration. 



from previous implementations o f the WSR algorithm which used a Landweber iteration with a more 
or less effective step size control ( Schriiver et al.J[2006l ). For these codes, the number of iterations was 
found to increase with ri} . Note, however, that these authors measured the number of steps until the 
magnitude of (5B fell below a given limit rather than the convergence speed (i.e., the fractional decrease 
of L per iteration step). 

Of course, the precision of the analytic model, calculated on the respective grid positions increases 
with higher grid resolution (see circles in the diagrams in Figure [5]). For the n = 32 grid this precision 
is reached already after 100 iterations with const = 10, while for n = 64 the much higher analytical 
precision seems to be out of reach. 

It is well known that it depends among other details also on the initial iterate B*^"^ how closely 
a desired solution can be approached by Krylov-type iteration schemes (e.g., Landweber, steepest 
descent or conjugate gradients). During the iteration, the search directions ^B'*^ build up a subspace 
of the Hilbert space for B and discrepancies between the iter ate B*^') and the solution can at best be 
eliminated inside this subspace, usually called Krylov space ( Saad .1120031 ). Differences between B(°) 
and the solution which fall out of the Krylov space cannot be corrected. A good choice of the initial 
B^^) therefore does not only save computation time but may be a necessity to reach a solution at 
all, at least for ill-posed problems for which the Krylov space remains limited. In Figure [6] we show a 
diagram similar to the n = 64 case in Figure O but here the initial B^°) was chosen to be the potential 
field inside V for the given normal component boundaries but with the tangential field components on 
dV expected for the final force-free field. Note that again the Lorentz force and the divergence of the 
initial B'^'^^ do not vanish because the final, non-potential tangential field boundary values on dV have 
been enforced. 

The diagram shows that with this improved initial field, it is no problem to reach the error bounds 
of the analytical model if const < 5. Note, however, a larger value for the constant, i.e., a greater bias 
towards the Lorentz force in ([3]) during the iteration does not necessarily produce a smaller final error 
of the Lorentz force term. 

For the GR scheme, the iterates B*^*) are always numerically divergence free and at any step they 
satisfy the normal boundary conditions. We therefore cannot start this iteration from B^^^ =0 but take 
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Fig. 7 Mean square Lorentz force < [j*-'^ xB*-'' p > in the course of the GR {left) and the WSR (right) iteration. 
For the WSR iteration the constant was chosen to be 2. The curves are the result for the n = 32 and 64 grid 
as marked. For the latter we show results for an initial B'°^ = (064) and B'°^ = Bpot (P64) except for the 
boundary values. The dashed lines denote the respective residual Lorentz force for the discretized Low and Lou 
(1990) solution Borg. 
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Fig. 8 Same as Figure [T] but the maximum Lorentz force maxlj*-'-* xB*-*'! in V for the GR iteration (left) and 
the WSR scheme (right). 



the potential field as a starting point instead. Since the information about the tangential boundary 
field is stored in the boundary values of a, the divergence of B^"^ and the current exactly vanish 
in contrast to the WSR initial field. 

The GR scheme converges much more rapidly than the WSR algorithm. However, the GR iteration 
scheme does not guarantee a continuous decrease of a certain object function, such as L for the WSR 
case. Less than 10 iterations were needed to eliminate the residual Lorentz forces to the level of the 
discrete Low and Lou solution (sec Figure [7] and [8]) . However, once this level has been reached, the 
Lorentz forces could not be lowered any further but even increased slightly again. This holds both for 
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Fig. 9 Mean square difference < (B''-* — Borg)^ > of the iterated field to the original Low and Lou (1990) 
solution for the GR iteration (left) and the WSR scheme (right). The coding of the curves refer to the same 
calculations as in Figures [7] and |S] 
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Fig. 10 Same as Figure [Q] but the maximum difference max|B''' — Borg| of the iterated field to the original 
Low and Lou (1990) solution for the GR iteration (left) and the WSR scheme {right). 



the mean and the maximum residual Lorentz force. Note that in Figs. [7] we show the mean square 
Lorentz force without weight w l/|i?pot| as in Figure[51 

From our calculations we may conclude therefore that with a WSR code and a potential field as 
initial iterate smaller residual Lorentz forces may eventually be reached than with the GR approach, 
provided the WSR code is supplied with consistent and exact boundary conditions of, e.g., a known 
analytic force free field. In this case a minimum L = exists and we are confident that our code will 
approach towards it due to its exact line search capabilities. The number of iterations though may be 
prohibitive in some cases. 

We have seen in the field line plots, that the GR result, even though it retains residual Lorentz forces, 
approaches the discrete Low and Lou solution extremely well. This can be confirmed if the difference 
B(«) — Borg is monitored during the iteration. As is shown in Figs. [H and [TU] the GR code here performs 
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much better than WSR. The mean square error it achieves is almost an order of magnitude smaUer than 
the WSR code can reach, even with the potential initial field. It seems that the most consistent field 
with the least Lorentz forces is reached after only 5 (n = 32) or ~ 7 (n = 64) iterations. Continuation 
of the iteration brings the field still closer to the Low and Lou solution but increases the level of the 
Lorentz force slightly. After about 10-14 iterations the difference B^*) — Borg is minimized and rises 
thereafter. A critical point for a practical application of the GR code is therefore when to stop the 
iteration. The saturation of the residual Lorentz forces seems to be a helpful indication. 

To measure differences of our final iterate with the original Low and Lou field, we calculate the 
following two error norms: 

S J2 \Bi - (Borg)i 

^ elements i£{x,y,z} 

J2 J2 l(Borg)i| 

elements i^{x.y ,z} 

J2 J2 \B^-{Borg)t\ 

elements i^{x,y,z} 

^""^ E E (Borg)f ■ 

elements i£{x,y,z} 

Ejj measures a normalized difference to the original field, Ec a correlation with it. The two error 
measures are not independent: for small errors, 1 - Ec ^ El,. Note that we do not subtract or 
correlate vectors but every individual component because from our grid, we obtain the values for 
different components at different locations. For the final GR result on a n = 64 grid, we obtain 

Ed = 3.81 10^3 ; Ec = 1.000034 

after 14 iteration steps which took 18 minutes of calculation on a 667 MHz Pentium IV Linux operated 
single processor with 256 MB Ram. The residual Lorentz forces of the reconstructed field had the same 
level as the discretized original Low and Lou field and the divergence was of the order of the numerical 
roundoff error. 

The equivalent WSR computations starting from B'"^ = Bpot yield 

Ed = 1.81 10-^ ; Ec = 0.99947 

after 640 iteration steps and 22 minutes of calculation time on the same computer. The mean square 
residual Lorentz forces of the reconstructed field could be reduced to about a third of those of the 
discretized original Low and Lou field, its divergence to 1/80 of < (V-Borg)^ >■ 



5.2 Twisted loop 

As a second experiment we tried to model a twisted magnetic loop with the a amplitude as a free 
parameter. The boundary conditions for the GR iterations were chosen in the following way: 

n-B = and a = on dV/{x \ z = 0} , 

v-^ (x — Xo-)^ 

z-B = yj crexp — and a = ainax|z-B| on dVr]{x\ z — 0} , (8) 

(T = ± 

where x^- = (0.5 — 0.2(7,0,0)-^ and r = 0.1. Since |z-B| has a maximum of 1, ttmax measures the 
maximum magnitude of a of our boundary values. We also rely on the symmetry of the boundary 
conditions and do not make the distinction between (dV)^ as is formally necessary (see eq.[2]). Rather, 
we would like to see our weighting mechanism (see section U) at work. All weights were set to unity. 
The grid size was chosen as n = 64. 

For the equivalent WSR calculations, we need to transform ^ into the full field vector on the 
domain boundary. However only the in-plane curl, n x V x B on dV is determined by the a boundary 
values. We therefore use the boundary fields from the GR results as boundary conditions for the WSR 
calculations. The local maximum discrepancy between a as prescribed in ([8]) and the GR result was 
less than 10%. 
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Fig. 11 Representative field lines of the reconstructions of the twisted loop model by the GR (left) and WSR 
(right) algorithms. In both cases Omax = 5 was chosen. The colour code at the bottom represents B^, the top 
plane shows the vertical projection of the field lines. The starting points of the field lines were chosen identically 
for both model results. 
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Fig. 12 Mean square Lorentz force (left) and magnetic field correction step (5B (right) for the twisted loop 
model reconstruction by the WSR code. The boundary conditions chosen correspond to Qmax ~ 5. 



In Figure [TT] we present the results of the different algorithms for a value ctmax = 5. A visual 
inspection reveals that the low lying loops correspond well in shape, however the outer loops from the 
GR reconstruction show much more twist than those produced with the WSR code. The reasons for 
this difference may be twofold: The WSR result has the larger residual Lorentz forces (see Figure [T^ 
and it took almost 2000 iterations until a level was reached from which < (j x B)^ > could not be 
lowered any further. In the course of this WSR iteration we observed that the twist of the outer field 
lines developed only very late. A continuation of the iteration may therefore bring a slight improvement 
towards the shape of the GR produced field lines. 
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The GR iterations, on the other hand, probably overtwist the outer loops. The reason is the 
integration of a along field lines in ([5|). On the ground surface the vertical current density aBz is 
confined to only two small concentrated spots ~ exp — 2((5x/r)^ of less than 10 grid spacings in diameter. 
Even though we take great care to follow as closely as possible the characteristics when we solve (O, 
a slight diffusion of a off the characteristics can probably not be avoided. The diffused current density 
probably causes a stronger twist for the outer field fines than is reafistic. 

The twisted loop model was designed to see how the codes behave when a is enhanced. We pursued 
this question only for the GR code. The amount of twi st with which a field line can be charged, i.e., 
the magnitude of amax is limited by the virial theorem ( Molodensky J[T969, ) . The virial theorem states 
that for a sphere of radius R 

J |B|2 ^rJ (|n-B|2 - InxBp) . (9) 

V av 

For a half space (i.e., R oo) this reduces to 

J (|i-B|2- |zxB|2) i 0. (10) 

z=0 

This is in fact one of th e relations enforced by the preprocessing scheme for the boundary data by 
Wiegelmann et al. (|2005l ). We can express the surface field in terms of a potential and a stream function 



z X B = z X V(p — V?/" , 

and it is clear that both cj) and ip will depend on the shape and amplitude of the a boundary condition 
at z = 0. If the boundary conditions are chosen as in then 

Ail;{a) = z-j = a(z-B) = a,„ax|z-BI(z-B) . 

The last step is due to the rigid connection between a and in ([5]). Hence the stream function 
is directly proportional to ckmax and invariable in its shape. To keep the balance in (|10|) . | V(/)(aniax)| 
has to decrease as |VV'(Q;max)l grows with amax- Setting V'(ci^max) = amaxV'(l) and making use of the 
orthogonality (due to the symmetry in ([8]) about x = 0.5) of the two fields z x V0 and 'Vtp on the 
plane z = we obtain from (jTU]): 

/ (|i?.|2-|V</>(amax)|2) / W 

2 ^ z=0 ^ z=0 



/ |V^(1)P - / |V^(1)P 

2=0 z=0 



A definite upper bound for lofmaxl is therefore reached when | V(/)(aniax)| has declined to 0. Numerically, 
we obtain for (g]): /^^^ = 0.031 and /^^g |VV'(1)P = 0.00011 which gives la^axl < 16 as upper 

bound. Since | V0(aniax)P probably never becomes zero, a realistic upper bound is much less than this 
crude estimate. 

In the GR iterations of the twisted loop model we varied |amax| in the range [0,10]. We observe that 
the convergence drastically decreases and comes to a halt between |amax| — 7 and 8 (see Figure [T5|) . 
We presume that this loss of convergence is not a failure of the code but it rather indicates the upper 
limit of feasible values of |amax|- Beyond this limit, a stationary force-free solution with a boundary as 
in ([8]) probably is not possible. Note that this transition is not so well visible in the residual Lorentz 
forces which continually rise as |aniax| and hence the general level of the current density increases. 

A point worthy of investigation is how the above estimate of the maximum current density of a 
stationary force free flux tube complies with the kink i nstability threshold for twisted flux tubes (e.g., 
iMikic. Schnack. and van Hovenj[l990l : IVan HovenllMokf ) . This instability criterium predicts stationarity 
only if the number of turns a fleld makes around the flux tube axis is less than about 2.4. 
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Fig. 13 Mean square Lorentz force (left) and magnetic field correction step 5B {right) for the twisted loop 
model with different amplitudes Qmax in the course of the GR iteration. T he respective values of Omax are 
attached to the curve. Note the different ordinate scales compared to Figure [121 



6 Discussion 



Force-free extrapolation c odes available nowada ys are capable to compute a field model in boxes as 
big as 256x256x256 (e.g.. ISchriiver et al.j|2006l ). Compared to the resolution of modern vector mag- 
netograms with pixel arrays of 1000x1000, the calculated extrapolation models are most often quite 
limited in their resolution. Keeping in mind that they have to oversample the observation in order to 
limit the discretization error, the obser ved data typically has to be smoothed bef ore it can serve as 
boundary condition for an extrapolation ( Wiegelmann. Inhester. and Sakurai.ll2005r ). There is therefore 
a definite need for fast and effective force-free reconstruction codes which can handle bigger models in 
order to make proper use of the resolution which the observations provide. 

In this paper we have implemented and tested two alternative schemes for the extrapolation of a 
force-free field from boundary data. The algorithms differ considerably and in order to compare these 
different approaches unobscured by differences in the numerical coding, we implemented them as far as 
possible in a similar way. The discretization on a finite element grid by means of Whitney forms results 
in codes with a very efficient performance. In addition, we have improved the WSR code considerably 
using an efficient conjugate gradient iteration. 

The two schemes differ not only in their approach towards a solution, but they also differ in the 
boundary information which has to be supplied. The WSR code requires boundary information which 
clearly overdetermines the problem unless parts of the boundary fields are left to be varied. But the 
full magnetic field vector on even part of the boundary is a very strong constraint and there is a great 
danger to impose inconsistent boundary values. This may be one reason w hy the WSR code converges 
well for known solutions like the Low and Lou model ( Low and LouJll990D where precise and consistent 
boundary values are supplied, while a little less reliable boundary values like those we retrieved from 
the result of a GR extrapolation slow down the WSR convergence speed markedly. 

The WSR code has as a free parameter the weight between the Lorentz force and the divergence 
term in ([3]). In our implementation it is hidden in the constant of the weight function w. We found const 
~ 1 a good value for the problems which we have delt with here. However, since the ([T]) is nonlinear, 
the performance and also the optimal parameters may vary with the problem studied. As an example, 
note the difference in the iteration history for the Low and Lou model for const — 2 and 5 in Figure [6l 
Here, the convergence changes drastically if the constant is modified by only a small amount. 
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Typical computation times on a n = 64 grid with our code are about 20 minutes on an ordinary 
home compiitcr. At present, the codes still include some checkout overhead, which could be dispensed 
with. Major improvements can probably be obtained by spreading the calculation onto multiple grids. 
We have made first tests with the WSR code for an adaptive enhancement of the grid resolution. Here, 
the solution on the lower grid n was used as initial iterate on the next finer grid 2n. We estimate that 
this way the computation time can be reduced by about a factor 1/3. Even more can probably be 
gained if the a true multigrid scheme is used. 
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